Fast numerical test of hyperbolic chaos 
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The effective numerical method is developed performing the test of the hyperbolicity of chaotic 
dynamics. The method employs ideas of algorithms for covariant Lyapunov vectors but avoids their 
explicit computation. The outcome is a distribution of a characteristic value which is bounded 
within the unit interval and whose zero indicate the presence of tangency between expanding and 
contracting subspaces. To perform the test one needs to solve several copies of equations for infinites- 
imal perturbations whose amount is equal to the sum of numbers of positive and zero Lyapunov 
exponents. Since for high-dimensional system this amount is normally much less then the full phase 
space dimension, this method provide the fast and memory saving way for numerical hyperbolicity 
test of such systems. 
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Introduction. - The dynamics is said to be uniformly 
hyperbolic if exponential rates of tangent space growth 
and contractions are always bounded and differ from zero 
by some global constants. Systems with hyperbolic dy- 
namics admit rigorous proof of their chaotic properties. 
The hyperbolic chaos is structurally stable, i.e., it per- 
sists under change of parameters of a system. Over 
many years this paradigm remained mainly the theoret- 
ical, since the most of the known systems do not con- 
form with the basic assumptions of uniform hyperbol- 
icity. But recently the interest to hyperbolic chaos has 
been renewed after a series of publications by Kuznetsov 
and his collaborators who suggested sufficiently simple 
ideas of practical implementation of the uniform hyper- 
bolic chaos in natural systems p], Hj . 

Tangent space at each point of a hyperbolic chaotic 
trajectory is split into invariant subspaces with differ- 
ent expanding and contracting properties. For discrete 
time systems there are two subspaces. The first one con- 
tains all expanding directions associated with positive 
Lyapunov exponents, and the second one consists of the 
all contracting directions corresponding to negative ex- 
ponents. The hyperbolicity implies the strict separation 
of these subspaces [3j . It means that the smallest angle 
between expanding and contracting vectors are globally 
bound from zero. For flows one more invariant neutral 
subspace is associated with zero Lyapunov exponents. To 
extend the definition of the hyperbolic dynamics to this 
case one have to require the strict separation of these 
three subspaces [4(. 

Rigorous mathematical verification of hyperbolicity is 
not always possible. To test this property numerically 
one can find bases for expanding, contracting (and neu- 
tral, if any) subspaces and compute the smallest angles 
between vectors from these subspaces for different points 
along trajectories. If the dynamics is non-hyperbolic, 
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among trajectories there are ones with tangencies, i.e., 
with zero angles between these vectors. Performing nu- 
merical simulations one normally can not hit such tra- 
jectory exactly. But one can expect that randomly cho- 
sen trajectory will pass infinitely close to the trajectories 
with tangencies. Thus for non-hyperbolic dynamics the 
distributions of smallest angles between subspaces have 
to be infinitely close to the origin, while for hyperbolic 
dynamics these distributions are well detached from zero. 

This approach was initially developed for low- 
dimensional dynamics (see [l[ for review) . Its application 
for high-dimensional dynamics became possible only re- 
cently after discovery of effective algorithms for compu- 
tation of covariant Lyapunov vectors (CLVs) [HQ. CLVs 
are associated with Lyapunov exponents and are covari- 
ant with the tangent flow, i.e., the z-th covariant vector 
at time t\ is mapped to the z-th covariant vector at t2- 
These vectors provide a natural bases for expanding, con- 
tracting, and neutral subspaces. 

In paper Q the angles between expanding and con- 
tracting subspaces spanned by CLVs were analysed to 
detect hyperbolic and non-hyperbolic regimes of a spa- 
tially extended system. Also the angles between sub- 
spaces spanned by CLVs were employed in Refs. [H, || to 
identify possible bases for inertial manifolds (lp| . 

Unfortunately, analysis of angles between subspaces 
spanned by CLVs is very resource-consuming. One needs 
to process a lot of angles for sufficiently long trajectory 
to obtain a representative statistics. Moreover, it is usu- 
ally unclear which CLVs from the negative end of the 
spectrum can be safely omitted without distortion of the 
result, and the most reliable way is to compute them all. 

The other quantitative numerical test of hyperbolic- 
ity is based on the direct verification of the cone cri- 
terion [J, This method is well developed for low- 
dimensional system and provides perhaps the most re- 
liable result. But unfortunately it is unclear yet how to 
transfer it to high-dimensional systems. Moreover even if 
this can be done, most probably it will also be a resource- 
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consuming. 

The paper reviews the theory of Lyapunov vec- 
tors and the related numerical algorithms. Moreover an 
improved approach for computing CLVs is suggested. In 
this letter the idea of this approach is employed to de- 
velop the fast method for testing of the hyperbolicity 
without explicit computing the CLVs. Given the sys- 
tem with fc_|- positive and ko zero Lyapunov exponents, 
it is enough to solve only k+ + ko equations for infinites- 
imal perturbations. Since for high-dimensional systems 
k + + k is normally much less then the whole phase space 
dimension, this method provides the fast and memory 
saving way for numerical hyperbolicity test of such sys- 
tems. 

Theory. Preforming the standard procedure of 

computation of Lyapunov exponents, we deal with a 
set of orthonormal vectors which are evolved in tangent 
space by alternating time evolution over an interval Tqr 
and re-orthogonalization via Gram-Schmidt or QR proce- 
dure. After sufficiently many iterations these vectors con- 
verge to backward Lyapunov vectors (also they are called 
Gram-Schmidt vectors). They are "backward" because 
initialized in the far past. Let ip~ (t) be the i-th backward 

Lyapunov vector and = [<pi (t), ¥>2 C0> • ■ • > VmO)] 

be an orthogonal matrix of these vectors, where m is 
the dimension of the tangent space. Performing the 
algorithm for Lyapunov exponents backward in time 
we obtain in the same way another set of orthonormal 
vectors <P\(t) and the corresponding orthogonal matrix 
*+ (i) = [cp+ (t), ip+ (f ), . . , , ¥>+(*)] . These vectors are re- 
ferred to as forward Lyapunov vectors, because they are 
initialized in the far future. Both f~(t) and <pf(t) are 
associated with the i-th Lyapunov exponent Aj. The Lya- 
punov exponents are assumed to be ordered as usual in 
descending order. 

Let r(t) = ["fi(t), 72(^)1 •■• ,7 m ] be a matrix of the 
CLVs 7j(i). These vectors are related to the forward and 
backward Lyapunov vectors as @, [Tl| 

r(t) = *-(t)A-(t) = S + (t)A+(f), (1) 

where A - (t) is an upper triangular matrix and A + (t) 
is a lower triangular matrix. In fact <I? - (t)A~(i) is the 
QR decomposition of T(t) while (t)A + (t) is its QL 
decomposition. 

Given the matrices <&~(t) and & + (t) one can com- 
pute covariant Lyapunov vectors using the equation 
P(*)A"(i) = A+(t), where 

P(t) = [*+(t)] T *-(t) (2) 

is a m x m orthogonal matrix, and A ± (t) are components 
of the LU decomposition of P(t). The details of this LU- 
method of computation of CLVs can be found in 

Let us assume that we have the whole set of m CLVs. 
Consider two subspaces of the tangent space. The first 
one is spanned by the first k CLVs, and the second one 
is spanned by the rest of them. We need to check if 
there is a tangency between vectors from these subspaces. 



Consider two arbitrary unit vectors from these subspaces: 

k m 

»i(*) = l> 7i (t), = E Wi®' ( 3 ) 

i=l i=k+l 

where and are some expansion coefficients. To detect 
the tangency we can compute the angle between V\{t) 
and i)2(t) and find such coefficients and C4 that mini- 
mize is. Zero minimal angle indicates the tangency 

As follows from Eq. fl} , the covariant vector 7 • (t) be- 
longs to the subspace spanned by the first j forward Lya- 
punov vectors cp~(t), where i — l,2,...,j. Also 7„(t) 
belongs to the subspace of the rest m — n + 1 forward 
vectors (p^(t), where i = n, n + 1, . . . , m. It means that 
V\{t) can be represented as a linear combination of the 
backward Lyapunov vectors and v^it) is a linear combi- 
nation of the forward Lyapunov vectors, 

k TCI 

v 1 (t)=J24<P7(^ «*(*)= E fat®, ( 4 ) 

i=l i=fe+l 

with some expansion coefficients e\ and c[. 

The tangency occurs if there are such e[ and that 
v 1 (t)=v 2 {t): 

k m 

E e ^ r (*)= E fat®- (5) 

i=l i=fe+l 

To find e\ and d i we multiply §5§ by vectors <Pj(t). The 
first k multiplications result in the set of equations for 
e\ and the rest of them produces the equations for 
provided that the nontrivial solution for e\ exists. The 
equations for e\ read: 

k 

£e^t(t) vr (t)=0, ./ 1.2..../,-. (6) 

i=l 

The nontrivial solution exists, when the scalar products 
<p^ (t)ip~ (t) , where 1 < i,j < k, form a singular matrix. 
One can see from Eq. ^ , that this matrix is the top left 
k x k submatrix of P. This provides an idea of the test 
of tangencies between subspaces and, in particular, the 
test of the hyperbolicity. 

Formulation of the method. - Given a dynamical sys- 
tem, in its tangent space we need to find the smallest 
angle between vectors from the subspace spanned by the 
first k covariant vectors and from the subspace spanned 
by the rest m — k of them, (a) We begin to move forward 
in time solving simultaneously the basic equations and 
k sets of linearized equations for infinitesimal perturba- 
tions. In the same way as in computing Lyapunov expo- 
nents, we alternate time evolution over intervals Tqr and 
orthogonalizations via QR or Gram-Schmidt algorithms. 
These steps are repeated for a while until columns of or- 
thogonalized matrices converge to tp~ (t) . (b) The same 
procedure is continued but now <p~(t) are saved. These 
steps are repeated as many times as many points of the 
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trajectory we are going to test, (c) We proceed to move 
forward along the trajectory with the basic system only 
and save the trajectory points. The equations for per- 
turbations are not solved on this stage, (d) We turn 
back and start to move backward along the saved tra- 
jectory performing steps with the basic system and k 
copies of equations for perturbations. Again time evo- 
lution alternates with orthogonalizations. In this way 
orthogonal vectors converge to <pf(t). There two subtle 
points here. First, the modified perturbation equations 
have to be used. The Jacobian matrix that determines 
these equations have to be transposed and its elements 
have to change their signs. In this case the vectors <pf(t) 
are computed in the proper order, i.e., k perturbation 
equations produce the k first vectors, and not the last 
ones (see [Il| for explanations). Second, dealing with 
the dissipative systems, it is better to use saved trajec- 
tory points avoiding solving the basic system backward 
in time. This is because negative Lyapunov exponents 
have typically large absolute values, and this can result 
in the dramatic instability of the numerical procedure, 
(e) After the arrival at the time step for which <f~(t) 
was previously saved, we proceed the above procedure, 
but additionally construct submatrices P(l :k, 1 : k), see 
Eq. and compute their normalized determinants as 

d fc = |det[P(l:M:fc)]|/fc! (?) 

Here < dk < 1 since absolute values of elements of P 
are less or equal to 1. (f) We accumulate many dk and 
compute their distribution. If this distribution is well 
separated from the origin, the trajectory never pass close 
to points of tangencies between two studied subspaces. If 
the amount of the processed points is sufficiently large, 
this is the reason to conjecture that the tangencies are 
absent at all. Otherwise, the close approach of the dis- 
tribution to the origin is the reason to conjecture that 
there are trajectories with the exact tangencies. 

Testing the hyperbolicity of chaos, we need to know 
the amount of positive k + and negative ko Lyapunov ex- 
ponents. If zero exponents are absent, the distribution of 
dk + have to be analysed. In presence of zero exponents, 
there are two distributions to consider: for dk + and for 
dk + +k - For hyperbolic chaotic systems both of the dis- 
tributions are well separated from the origin. Since the 
number of expanding and neutral directions are normally 
much less then the full dimension of the tangent space, 
the represented algorithm provide sufficiently fast way of 
testing the hyperbolicity even for high-dimensional sys- 
tems. 

Examples. - First we consider the Lorenz system 

x = a(y — x),y = rx — xz — y,z = xy — bz. (8) 

where r = 28, b — 8/3, a = 10. The Lyapunov exponents 
are 0.90, 0, and —14.57. There is one positive and one 
zero exponent and thus we have to analyze distributions 
of g?i and g?2, see Fig. [TJ This system is known to be sin- 
gular hyperbolic, which means its tangent space admits 
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Figure 1. (color online). Distributions of normalized de- 
terminants (0 over a trajectory of the Lorenz system (J8|. 
T QR = 0.1, At « 0.002. 

an invariant splitting E c © E ne into a 1-dimensional uni- 
formly contracting sub-space and 2-dimensional volume- 
expanding subspace [12]. In an agreement with this defi- 
nition in Fig.[T]the distribution p(d±) touches the origin, 
while p(da) does not. It indicates the existence of trajec- 
tories with tangencies between the expanding and neu- 
tral subspaces, while the contracting subspace remains 
isolated. 

The other example is the system of two alternately 
excited van der Pol oscillators 

x — [Acos(27rf/T) — x 2 ]x + ujqX — eycosuiot, 
ij - [-Acos{2iTt/T) - y 2 ]y + 4cu 2 y = ex 2 . 

where A ~ 5, T — 6, e — 0.5, ojq — 2tt. Corresponding to 
this system stroboscopic map at t = t n = nT is known 
to be hyperbolic P, 0] • The Lyapunov exponents for this 
map are 0.68, —2.61, —4.61, —6.10. There is one positive 
exponent and the zero one is absent since the system is 
non-autonomous. Hence the hyperbolicity test for this 
map includes the analysis of p{d±). As we see in Fig. [2] 
values of d\ are located far from the origin, confirming 
the applicability of the method. 

Next we consider the complex Ginzburg-Landau equa- 
tion 

d t a = a- (\ + ic)\a\ 2 a+{l + \b)d 2 x a (10) 

at c = 3, b = —2, that corresponds to the amplitude 
chaos [I3j. Since the no-flux boundary conditions are 
used, the system has two continues symmetries, i.e., time 
translation and phase rotation, and thus has two zero 
Lyapunov exponents [Til ]. For chosen parameter values 
the first seven exponents read: 0.30, 0.18, 0.060, 0, 0, 
—0.085, —0.27. It means that we have to consider p{d%) 
and p(ds), see Fig. [3] For both distributions dk are es- 
sentially nonzero at the origin, thus we have to conclude 
that the chaos is essentially non-hyperbolic. 

Finally, we apply the test to the alternately excited 
Ginzburg-Landau equations, which are the amplitude 
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Figure 2. (color online). Distribution of dk for the strobo- 
scopic map built for the system © at t = t n = nT. Tqr = T, 
At « 0.01. 
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where A = 3, T = 5, e = 0.05. The spatio-temporal 
chaos in this system is hyperbolic when its length L is 
sufficiently small and the hyperbolicity is destroyed when 
the length growths 0] • We consider the same parameters 
and numerical mesh as in Q. For L = 15 first five Lya- 
punov exponents, computed for the stroboscopic map at 
t = t„, = nT, are 0.69, 0.69, -0.11, -0.61, -1.62. For 
L = 17 the exponents are 0.69, 0.67, 0.13, -0.28, -0.85. 
Hence we have to consider p(cfe) for the first case and 
p{d^) for the second case, see Fig. 21 One can see that 
represented method detects the reported in Ref. @] tran- 
sition from the hyperbolic chaos at L — 15 to the non- 
hyperbolic one at L = 17. 

Altogether we see that the demonstrated examples 
agree well with the previously known results. The 
method works sufficiently fast. All distributions repre- 
sented above have been computed for 10 5 values of dk 
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Figure 3. (color online). Distributions of (a) d% and (b) ds 
for complex Ginzburg-Landau equation Tqr = 0.5, N = 
128, Aa; = 0.1, At w 0.0005. 



Figure 4. (color online). Distribution of for alternately 
excited Ginzburg-Landau equation (|11[1 . Tqr = T, N — 51, 
Ax = L/(N — 1), At w 0.05. 



equations for the system © supplied with diffusive usmg 2.6 GHz processor. For distributed systems it took 

terms: approximately 3 hours. 

8 t a = Acos(27rf/T)a - \a\ 2 a - ieb + dla, The research is supported by RFBR-DFG grant 11-02- 

(11) 91334 

dtb = -Acos(2nt/T)b - |6| 2 & - iea 2 + d 2 x b, 
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